An Improved Gauge Driver for the Generalized Harmonic Einstein System 
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A new gauge driver is introduced for the generalized harmonic (GH) representation of Einstein's 
equation. This new driver allows a rather general class of gauge conditions to be implemented in a 
way that maintains the hyperbolicity of the combined evolution system. This driver is more stable 
and effective, and unlike previous drivers, allows stable evolutions using the dual-frame evolution 
technique. Appropriate boundary conditions for this new gauge driver are constructed, and a new 
boundary condition for the "gauge" components of the spacetime metric in the GH Einstein system 
is introduced. The stability and effectiveness of this new gauge driver are demonstrated through 
numerical tests, which impose a new damped-wave gauge condition on the evolutions of single 
black-hole spacetimes. 



I. INTRODUCTION 

The gauge (or coordinate) degrees of freedom in the 
generalized harmonic (GH) form of the Einstein equa- 
tions are determined by specifying the gauge-source func- 
tions iJ". These functions are defined as the results of 
the CO- variant scalar- wave operator acting on each of the 
spacetime coordinates x"". 



(1) 



(We use Latin letters from the beginning of the alpha- 
bet, a, b, c, for spacetime indices.) The GH form of 
Einstein's equations can be represented (somewhat ab- 
stractly) as 

^'^dcddA^b + daHf, + dbHa = Qab{H, V, d4>), (2) 

where ipab is the spacetime metric, Ha ~ ipabH^i and 
Qab represents lower-order terms that depend on Ha, 
the metric, and its first derivatives. These equations 
are manifestly hyperbolic whenever Ha is specified as 
an explicit function of the coordinates and the metric: 
Ha = Haixjtp)- In this case the terms daHi, appearing 
in Eq. ([2]) contain at most first derivatives of the met- 
ric. The Einstein equations become, therefore, a set of 
second-order wave equations for each component of the 
spacetime metric: 



tp" dcddlpab Qabix, Tp, dljj) . 



(3) 



Thus the Einstein equations are manifestly hyperbolic for 
any Ha = Ha{x,'p). 

Most of the useful gauge conditions developed by the 
numerical relativity community over the past several 
decades can not, unfortunately, be expressed in the sim- 
ple form Ha = Ha{x, ip), unless the full spacetime metric 
ipab ~ i^ab{x) is known a priori. Many of these condi- 
tions (e.g., Bona-Masso slicing or the T-driver shift con- 
ditions) would require gauge-source functions that de- 
pend on the spacetime metric and its first derivatives: 
Ha ~ H a{x , ip , dip) , cf. Ref. In this case the terms 
daHb in Eq. ([2]) would depend on the second derivatives 
of the metric, ipab, and this (generically) destroys the 
hyperbolicity of the system. 



This problem can be overcome by elevating Ha to the 
status of an independent dynamical field and introduc- 
ing suitable evolution equations for Ha, which we call 
gauge drivers [H, 0, Q- One obvious choice is to con- 
struct gauge-driver equations that force Ha to evolve to- 
ward the desired gauge, e.g.. Ha — > Fa where Fa is the 
target for the selected gauge. To be useful these gauge 
driver equations must also make the combined Einstein 
gauge-driver system hyperbolic. It is fairly easy to con- 
struct hyperbolic evolution systems designed to evolve 
Ha toward any target Fa{x,'ip,dip) that depends on the 
spacetime metric and its first derivatives Many of 
the gauge conditions found most useful by the numerical 
relativity community have targets Fa that belong to this 
class. In most cases however, the coupled Einstein gauge- 
driver evolution equations are unstable and the evolved 
Ha does not evolve robustly toward every target Fa in 
this class for generic evolutions. The Einstein gauge- 
driver system is very complicated, and there are many 
opportunities for unstable couplings to develop between 
the dynamics of the spacetime metric and the dynam- 
ics of the gauge field Ha- Some gauge conditions, in- 
cluding certain Bona-Masso slicing conditions and some 
versions of the F-driver shift conditions, have been imple- 
mented fairly successfully using gauge drivers of this type 
in full 3D evolutions of strongly perturbed single black- 
hole spacetimes However, we find that even these 
"successful" gauge drivers fail when more complicated 
simulations are attempted, e.g., evolving a single black 
hole in a rotating reference frame or evolving black-hole 
binary systems. 

The purpose of this paper is to develop a better gauge 
driver that overcomes some of these problems. To this 
end we introduce in Sec. |ll] a new class of "first-order" 
gauge driver evolution equations, which are considerably 
simpler than earlier drivers. The dynamical simplicity 
of these new drivers reduces the internal dynamical de- 
grees of freedom available to Ha (in a sense discussed 
in more detail in Sec. |ll|, hence reducing the possibility 
of unwanted feedback or resonances with the dynamics 
of the Einstein system. We describe numerical tests of 
this new gauge-driver system in Sec. IIIII that use a new 
damped-wave gauge introduced in Appendix [X] to pro- 
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vide an interesting non-trivial dynamical target Fa- Us- 
ing this target Fa we perform a series of numerical tests 
that evolve single black-hole spacetimes with large dy- 
namical gauge perturbations. These tests demonstrate 
the effectiveness and stability of the new gauge-driver 
system for single- and dual-coordinate frame evolutions. 
The strongly perturbed black holes in these tests always 
evolved into non-singular time independent states, which 
suggests that the new damped- wave gauge conditions in- 
troduced here may prove to be useful for numerical simu- 
lations of more general dynamical black-hole spacetimes 
as well. 

We describe in some detail a number of technical prop- 
erties of this new gauge-driver system in a series of Ap- 
pendices. In Appendix |B] we show that any member of 
this new class of first-order gauge drivers can be cou- 
pled to the GH Einstein system in a way that makes the 
combined system symmetric hyperbolic. In Appendix [Cl 
we develop a dual-coordinate frame version of this new 
gauge-driver system, which is needed to evolve black-hole 
binary systems for example. In Appendix [D] we analyze 
the evolution of the constraints in the new combined GH 
Einstein gauge-driver system. We show that the con- 
straints and their evolution equations are the same as 
those of the pure GH Einstein system, hence the con- 
straint damping properties of the original GH Einstein 
system are also unchanged. In Appendix lEl we construct 
boundary conditions for the gauge-driver system. In 
most cases these boundary conditions turn out to be the 
same as those used for the pure GH Einstein system, but 
their representations in terms of the characteristic fields 
of the gauge-driver system are different in some cases. 
We also introduce a new constraint-preserving boundary 
condition for the "gauge" components of the spacetimc 
metric in the GH Einstein system. 

II. FIRST-ORDER GAUGE DRIVER 

The gauge drivers previously introduced for the GH 
Einstein system [H [1, 01 were constructed by elevating 
the gauge-source function Ha to the status of a dynamical 
field that is evolved by a second-order wave equation for 
Ha having the general form, 

r^dcddHa = Qa{H, dH, V-, d^). (4) 

When this type of evolution equation for Ha is used to- 
gether with the GH Einstein evolution Eq. ([2]), the com- 
bined system is manifestly hyperbolic. The first imple- 
mentations of this type of gauge driver were fairly suc- 
cessful, allowing a few successful binary black-hole inspi- 
ral, merger and ringdown simulations 0, 0]- ^ disad- 
vantage of these first gauge drivers however is that they 
were not designed to drive Ha toward a predetermined 
target Fa, so using them made it difficult or impossible 
to predict what gauge would ultimately be imposed on 
the solution. One reason for this ambiguity is the dy- 
namical complexity of the operator used to evolve Ha- 



Even the homogeneous driver, Eq. ([4]) with Qa = 0, has 
a wealth of solutions that arc not naturally attracted to- 
ward any particular target Fa . So it is not surprising that 
these first gauge drivers have not been found to be very 
effective for implementing pre-determined gauge condi- 
tions or for performing evolutions in generic situations. 
The goal here is to introduce a gauge driver that drives 
Ha toward a predetermined gauge specified by Fa more 
robustly and in more generic situations than was possi- 
ble with the first gauge drivers of this type [l| based on 
the the complicated second-order wave operator used in 
Eq. Q. 

An ideal gauge-driver would determine Ha from an 
evolution equation like, 

dtHa^-KHa-Fa), (5) 

whose solutions all approach the target gauge-source 
function Fa exponentially, at a rate determined by the 
freely specifiable parameter fi. Unfortunately the evolu- 
tion system formed by combining Eq. ([5]) with the GH 
Einstein evolution Eq. ([2]), does not appear to be hyper- 
bolic. There is a simple generalization of this ideal gauge 
driver however that can be used with the GH Einstein 
equations to construct a composite evolution system that 
is hyperbolic. Let i° denote the future-directed normal to 
the constant-t hypersurfaces. Then the first-order gauge 
driver, 

t''dbHa = -KHa - Fa), (6) 

combined with the GH Einstein evolution Eq. ([2]) turns 
out to be a hyperbolic system. 

We present a proof below that the combined GH Ein- 
stein gauge-driver system, Eqs. ([2]) and is hyperbolic. 
Before turning to that technical issue in AppendixlBlhow- 
ever, we point out that the very simple gauge driver, 
Eq. ([6|), has some limitations which can be overcome 
to some extent by a simple modification. To see these 
limitations we introduce spacetime coordinates, {i, x*}, 
where the time coordinate t labels the leaves in a folia- 
tion of spacelike hypersurfaces on which the points are 
identified by the spatial coordinates x*. In this coordi- 
nate system we use the standard 3+1 representation of 
the spacetime metric, ipab- 

ds^ = tpabdx'^dx'', 

= -N^dt^ + g,j{dx' + N'dt){dx^ + N^dt), (7) 

where gij is the intrinsic spatial metric of the constant-i 
hypersurfaces, and N and A^* are referred to as the lapse 
and shift respectively. (We use Latin letters from the 
middle of the alphabet, i, j, k, for purely spatial in- 
dices.) The unit normal to the constant-i hypersurfaces, 
t", has the 3+1 representation fda = N-'^{dt - N''dk) 
in this notation. Thus the gauge driver given in Eq. ([6]) 
can be written more explicitly in 3+1 form as 

dtHa-N''dkHa = -fliHa-Fa), (8) 
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where fi = jlN . This gauge driver has the property that 
Ha is driven toward Fa as seen by observers moving along 
the world lines of the hypersurface normal t° . However at 
a fixed spatial coordinate, x', the quantity Ha — Fa is not 
necessarily driven to zero. Therefore the evolution of a 
dynamical spacetime (e.g., a perturbed black hole) using 
this driver will not evolve toward a time independent 
state in which Ha = Fa- Rather this driver will tend to 
evolve solutions into states with N^dkHa = fi{Ha ~ Fa). 
This gauge may provide a reasonable representation of 
the spacetime, but it will not be the gauge Ha = Fa the 
driver was intended to enforce. 

This limitation in the gauge driver of Eq. ([6]) can be 
overcome by introducing an additional dynamical field, 
da defined as 

dtea + rjea = -vN'^dkHa. (9) 
or equivalently, 

Oa{t) = -77 /* e'^^''-'^N''dkHa{t')dt'. (10) 

— oo 

The 6a field is an exponentially weighted time average of 
—N'^dkHa^hich can be used to modify the gauge driver 
of Eq. ® H: 

dtHa - N'^dkHa = -li{Ha - Fa) + 9a. (11) 

All time independent solutions of the first-order gauge 
driver consisting of Eqs. ^ and ([TT|) must now satisfy 
the desired gauge condition Ha = Fa. Since the gauge- 
driving parameters rj and fj, are freely specifiable, they 
can be chosen to enforce the desired gauge on a timescale 
shorter than the characteristic time r on which the space- 
time evolves. Thus we expect the desired gauge can 
be enforced using this driver with reasonable accuracy 
Ha ~ Fa in any spacetime. 

In Sec, mil we present numerical tests of this first-order 
gauge driver that demonstrate how well it succeeds. In a 
series of Appendices we also present some formal analy- 
ses of a variety of mathematical properties of the new 
gauge driver composed of Eqs. ^ and PT|) together 
with the GH Einstein Eq. In particular we show 

in Appendix |B] that this combined GH Einstein gauge- 
driver system is symmetric hyperbolic. In Appendix [Clwe 
construct a dual-coordinate frame version of this gauge 
driver that can be used for example in the evolution of 
binary black- hole spacctimes. In Appendix ID] we analyze 
the constraints and the evolution of the constraints in the 
GH Einstein gauge-driver system. And in Appendix lEl we 
formulate boundary conditions for the new gauge-driver 
system. 

III. NUMERICAL TESTS 



These tests evolve a Schwarzschild black hole with per- 
turbed lapse and shift using the full coupled non-linear 
equations for the GH Einstein gauge-driver system, as 
described in Sec. |TT1 We measure the stability and effec- 
tiveness of the new gauge-driver system as it attempts to 
drive this single black-hole spacetime from the isotropic 
maximal-slicing gauge used to specify the initial data to 
an interesting new damped- wave gauge introduced in Ap- 
pendix 13 

These numerical tests are conducted using the infra- 
structure of the Caltech/Cornell Spectral Einstein Code 
(SpEC). This code uses pseudo-spectral collocation 
methods, as described for example in Refs. d, Q. We 
use the generalized harmonic form of the Einstein equa- 
tions, as described in Ref. 0, together with the new 
gauge driver Eqs. ^ and pT|) . Some of the tests re- 
ported here use the dual-coordinate frame version of the 
new gauge-driver system described in Appendix [C] For 
these dual-frame tests we use the static Schwarzschild 
coordinates as the "inertial" frame, and a "co-moving" 
frame that rotates uniformly at angular velocity Q with 
respect to the inertial frame. The evolution equations for 
the combined GH Einstein gauge-driver system are inte- 
grated in time using the method of lines and the adaptive 
fifth-order Dormand-Prince integrator Q . 

Initial conditions are needed for any evolution of the 
combined GH Einstein gauge-driver system. These initial 
data consist of the spacetime metric ipab, its time deriva- 
tive dttpab, the gauge-source function Ha, and the time 
averaging field 9a. For the tests described here we take 
the initial spacetime metric ipab to be the Schwarzschild 
geometry plus perturbations as described below. Wc set 
the time derivatives of the spatial components of the met- 
ric initially to zero, and the time-derivatives of the lapse 
and shift, dtN and dtN^ , are chosen to make TV and TV' 
initially time independent. For the dual- frame evolution 
tests described below, these time derivatives are chosen 
to make N and the co-moving frame components of TV' 
time independent initially in the co-moving frame. The 
initial value of Ha is chosen to enforce the gauge con- 
straint, Ca = Ha+Ta = 0, initially. The value of the time 
averaging field 9a is set initially to ensure that its time 
derivative vanishes, as determined by Eq. ^ or (jC6j) . 

For these tests we construct initial data consisting of a 
Schwarzschild black hole with perturbations in the lapse 
and shift. For the unperturbed hole we use isotropic 
spatial coordinates and maximal time slices 0, [l0| . The 
unperturbed spatial metric in this representation is given 
by, 

ds^ = g^jdx'dx^ = {dx^ + dy^ + dz^) , (12) 

where — x^ + + , and R{r) (the areal radius) 
satisfies the differential equation. 



In this section we describe the results of 3D numer- 
ical tests of the new GH Einstein gauge-driver system. 
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The constant M is the mass of the hole, and C is a param- 
eter that specifies the particular maximal slicing. Finally, 
the unperturbed lapse N and shift iV* for this represen- 
tation of Schwarzschild are given by, 



N 



2M 



R 

Cf' / 2M 



(14) 
(15) 



where f* is the outward directed radial unit vector: 
gtjf'f^ = 1. 

We perturb this spacetimc by changing the initial val- 
ues of the lapse and shift, and their time derivatives. This 
type of perturbation changes the spacetime coordinates 
(or gauge) of the solution, but not its geometry. For these 
tests we modify the lapse and shift of Eqs. p4|) and (fTS]) 
by adding perturbations of the form, 

5N = Asin(2^r/ro)e-^'-'-<=)'/'"Vto„ (16) 
6N' = Asm{2^:r/ro)e-^''-'^'■-^'/'"\l„^f\ (17) 

where Yi„i is the standard scalar spherical harmonic. In 
our numerical tests wc use the background metric with 
C = 1.73il/^, and perturbations with A — 0.01, rc = 
15M, w = 3Af , ro = 6M, and ? = 2, m 0. 

These numerical tests are performed using the target 
gauge-source function for the new damped- wave gauge. 



Fa 



ML log ( ^ ) 



(18) 



where /i^ and us are damping parameters, g = det gtj, 
andp is a constant. This new gauge condition is discussed 
in some detail in Appendix |Al The gauge used to prepare 
the perturbed Schwarzschild initial data, Eqs. (fT^ - (fT7)) . 
is very different from the dampcd-wave gauge condition. 
It is always difficult to start evolutions in a smooth and 
convergent way using initial data prepared with a signif- 
icantly different gauge. To minimize this start-up prob- 
lem, it is common practice to turn on the new gauge 
condition gradually. We do this in our gauge driver sys- 
tem by defining an initial target Fi""* that is simply the 
constraint-satisfying Ha of the unperturbed initial data. 
Except for the perturbation, this is exactly the gauge 
needed for a time independent evolution of these initial 
data. We then set the target Fa to 



F, 



; = + (l - e-*'/^') i^i^^. 



(19) 



where F^"*^ is the target gauge-source function for the 
damped- wave gauge defined in Eq. p8|) . This choice 
for Fa changes the gauge condition from its initial state 
Fi"' to the desired F^^^ smoothly and gradually on the 
timescale T. For the tests discussed here we use T = lOM 
for the value of this time-blending parameter. 

These tests use the dampcd-wave gauge condition de- 
fined in Eq. (|18p with damping parameters /xg = fJ-L = 



0.1 and p = 0.5. Most of these tests (except as noted 
below) use the values ^ = 77 = 16 for the gauge-driver 
parameters, used in Eqs. Q and (|lip . and the boundary 
gauge-driver parameter = 1 used in Eq. (jElOp . These 
tests set the constraint damping parameters of the GH 
Einstein system to the values: 70 = 72 = 2 and 71 = — 1, 
cf. Rcf.fl. 

Wc perform these numerical tests on a computational 
domain consisting of a spherical shell that extends from 
r = 0.78M (just inside the horizon in the initial coor- 
dinates) to r = 60M (well outside the domain of influ- 
ence of the initial perturbations). We divide this domain 
into sixteen sub-domains, which allows us to distribute 
the computation over several processors to enhance com- 
putational speed. In each sub-domain we express each 
Cartesian component of each dynamical field as a sum of 
Chebyshev polynomials of r (through order TV,. — 1) mul- 
tiplied by scalar spherical harmonics (through order L). 
The radii of the inner and outer edges of the various sub- 
domains are adjusted to distribute the truncation error 
during the full time evolution more or less uniformly on 
the grid. The specific radii of the sub-domain boundaries 
used in these tests are 0.78Af, 1.68M, and k x 4.0M for 
k = l, ...,15. 

In the pseudo-spectral numerical method used here, 
each Cartesian component of each dynamical field is ex- 
panded as a sum of the form: 



1 L 



uir,e,^)= UMmn{r)Ye,^{e,^), (20) 

fc=0 e=o m=-l 

where the Ukem are referred to as the spectral coefficients 
of the field u. These spectral coefficients must be mod- 
ified in this method through a process called spectral 
filtering. We use two types of spectral filtering in these 
tests. One type affects the angular spectral coefficients, 
as described in Ref. This filter sets to zero in each 
time step the changes in the top four tensor spherical 
harmonic expansion coefficients of each of the dynam- 
ical fields. This filtering step is needed to eliminate an 
instability associated with the inconsistent mixing of ten- 
sor spherical harmonics whenever angular derivatives are 
computed in our approach. In addition we also perform 
the following radial filtering. 



1)]' 



Ukeri 



(21) 



where T{uktrn) represents the filtered coefficients, before 
applying outer boundary conditions as described in Ap- 
pendix [E] For these tests we use p = 0.9 and p = 18, 
which leaves essentially unchanged the coefficients Ukim 
with k < 2{Nr — l)/3, while the coefficient of the highest 
mode, k ~ Nr — 1, is effectively set to zero. This radial 
filter implements in a smooth way the standard 2/3 filter 
often used to cure non-linear aliasing that can occur in 
spectral evolutions [ll|, [l^l ■ 

The damped- wave gauge conditions defined by Eq. p8p 
(and described in Appendix [A]) are significantly different 
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FIG. 1: Coordinate radius of the apparent horizon of the black 
hole Rh as it evolves under the effects of the dynamically 
driven gauge. This test uses a single-frame evolution with 
gauge-driver parameters n = rj = 16. 



than those satisfied by the perturbed maximally sliced 
representation of the Schwarzschild geometry used as ini- 
tial data for this test. Consequently the representation 
of the black hole in our test becomes very dynamical, 
primarily due to these gauge differences, and also due to 
the presence of the asymmetric perturbation applied to 
the lapse and shift. Figure [1] illustrates just how signifi- 
cant these gauge differences are by showing the evolution 
of the coordinate radius of the apparent horizon Rh of 
the black hole. In these tests the radius of the apparent 
horizon Rh grows by 50%, changing from an initial value 
of 0.86M to a final radius of 1.28A/. 

Figure [5] illustrates the constraint violations for a 
single-frame evolution of the GH Einstein gauge-driver 
system, and demonstrates the stability and convergence 
of our numerical method. The constraints of the GH 
Einstein gauge-driver system are identical to those of the 
GH Einstein system, as discussed in some detail in Ap- 
pendix [d1 Therefore we measure constraint violations 
using the quantity || Cgh||, the ratio of an norm of all 
the GH Einstein constraints divided by an norm of 
the derivatives of the dynamical fields. This constraint 
norm vanishes iff the constraints are satisfied, and has 
been normalized to be of order unity when constraint vi- 
olations begin to dominate the solution. This constraint 
norm was originally introduced to measure constraint vi- 
olations for the pure GH Einstein system in Eq. (71) of 
Ref. The constraint violations become largest and 
the rate of convergence of the simulations decreases dur- 
ing the time interval 15M < t < 30M in Fig. [H when 
the inward moving gauge perturbation interacts most 
strongly with the black hole. These results show that 
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FIG. 2: Constraints of the GH Einstein system ||Cgh|| for 
a single-frame evolution of a Schwarzschild black hole with 
strongly perturbed lapse and shift. This test uses a single- 
frame evolution with gauge-driver parameters jj, = rj — 16, 
and several different values of the numerical resolution pa- 
rameters Nr and L. The small inset graph contains a magni- 
fied view of II CghII during the time interval 15M < t < 30M, 
showing that the solution is convergent during this most dy- 
namical part of the evolution. 



the constraints are well satisfied throughout the evolu- 
tions, demonstrates that our numerical methods are con- 
vergent, and shows that the GH Einstein gauge-driver 
system is stable over many dynamical timescales. Fig- 
ure [3] provides another illustration of the stability and the 
numerical convergence of the GH Einstein gauge-driver 
system. In this figure we show \dM{t)\/M the evolution 
of the difference between the evolved and the initial mass 
of the black hole (as determined from the area of its ap- 
parent horizon). 

Figure |4] demonstrates the effectiveness of the gauge- 
driver system for this test problem. The difference be- 
tween the gauge source function Ha and the target func- 
tion to which it is being driven. Fa, is measured using 
the following L" norm: 



\H - F\ 



j^m''\Ha-Fa){Hb-Fb)d'' 



j^m'^'^F^Fdd- 



X 



(22) 



where ra""^ is a positive definite matrix, set to the identity, 
rn"'' = 5°-^ , for these tests. This norm vanishes if and only 
if the target gauge condition. Ha = Fa is satisfied, and it 
is scaled so that Ha bears little resemblance to the target 
Fa whenever it becomes of order unity. Figure [5] shows 
that the initial mismatch between the gauge of the per- 
turbed black hole and the damped wave gauge conditions 
(defined by Fa) causes — i^||/||F|| to grow initially. 
But the gauge-driver steps in and limits this growth to 
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FIG. 3: Curves show |(5M|/M, the deviations in the mass of 
the hole from its initial value. This test uses a single-frame 
evolution with gauge-driver parameters /i = = 16. 
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FIG. 4: Effectiveness of the gauge-driver equation is demon- 
strated by showing \\H — F\\I\\F\\ for an evolution of a 
Schwarzschild black hole with strongly perturbed lapse and 
shift. This test uses a single-frame evolution with gauge- 
driver parameters = 77 = 16. 



a maximum of about 0.02 in these evolutions, and then 
drives — J^||/||i^|| to very small values (depending on 
the numerical resolution) at late times. 

The evolution tests illustrated in Figs. [TJ- [4] were per- 
formed using the single-frame version of the gauge-driver 
system described in Sec. Binary black-hole simula- 
tions are done with the Caltech/CorncU SpEC code using 



a dual-coordinate frame formulation of the GH Einstein 
equations . In this formulation the components of the 
various tensor fields are defined with respect to a non- 
rotating inertial coordinate frame, while the equations 
for these field components are solved using a co-moving 
coordinate frame that tracks the motions of the black 
holes. A dual-frame version of the GH Einstein gauge- 
driver system is developed in Appendix [Cl We have per- 
formed the same perturbed single black-hole evolution 
tests illustrated in Figs.[T]-[31 using this dual- frame ver- 
sion of the GH Einstein gauge-driver system. For these 
tests we use a co-moving frame that rotates with re- 
spect to the asymptotic inertial frame at angular velocity 

= XjM. (This means that equatorial grid points in this 
test move at 60 times the speed of light at the outer edge 
of our computational domain.) The gauge driver used 
for these evolutions is the hybrid driver described in Ap- 
pendix O Eqs. (jClOp and (|C11|) . This driver attempts 
to enforce the comoving-frame gauge condition iJ^ = F„^ 
in the spacetime region near the black hole, while en- 
forcing the inertial-frame condition Ha = near the 
outer boundary of the computational domain. The tran- 
sition between these is accomplished by smoothly blend- 
ing the two conditions at intermediate points using a 
weight function w(x), cf. Eqs. (|C10[) and (|Glip . In re- 
gions where w{x) = 1, the pure comoving-frame condi- 
tion is enforced, and where w{x) = the pure inertial- 
frame condition is used. For these numerical tests we use 
w{r) = f.-[r/(o.mRa)]" ^ ^here Ro = 60M is the outer 
radius of the computational domain. This choice accu- 
rately enforces the comoving-frame condition in the inner 
region of the domain where r < 2Ro/3, and the inertial- 
frame condition at points located very near the outer 
boundary, r ^ Rq. 

The graphs of the quantities depicted in Figs. [T]- [4] for 
the dual-frame evolution case are almost identical to their 
single-frame evolution counterparts. So we will not show 
those graphs again here. Instead we show in Fig. [5] a se- 
ries of evolutions performed with the dual-frame system 
in which the effects of varying the gauge-driver param- 
eters ^ and 77 are examined. We see from these results, 
that the gauge-driver system is very effective in driving 
Ha Fa for a wide range of gauge-driver parameters. 
Evolutions using larger values of the gauge-driver param- 
eters are generally more effective in keeping the quantity 
— F||/||F|| small and driving it quickly toward zero. 
The gauge-driver system is stable and effective over a 
rather wide range of parameters, but becomes ineffec- 
tive when the gauge-driver parameters get smaller than 
about one, and the system also becomes unstable when 
the parameters are larger than a few hundred. 



APPENDIX A: DAMPED- WAVE GAUGE 
CONDITIONS 

Harmonic gauge is defined by the condition that each 
coordinate x° satisfies the co-variant scalar wave equa- 
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FIG. 5: Effectiveness of the gauge-driver system is demon- 
strated for various values of the gauge-driver parameters rj 
and jj.. This test uses a dual- frame evolution method with the 
co-moving frame rotating with respect to the inertial frame 
at angular velocity Q — 1/M. 



tion: 



= = 0. 



(Al) 



Harmonic coordinates have proven to be extremely use- 
ful for analytical studies of the Einstein equations, but 
have found only limited success in numerical problems 
like simulations of complicated highly dynamical black- 
hole mergers. A likely reason for some of these difficul- 
ties is the wealth of "interesting" dynamical solutions 
to the harmonic gauge condition itself. Eq. (jAip . Since 
all "physical" dynamical fields arc expressed in terms of 
the coordinates, an ideal gauge condition would limit co- 
ordinates to those that are simple, straightforward, de- 
pendable, and non-singular; having "interesting" dynam- 
ics of their own is not a desirable feature for coordinates. 
We propose to reduce the dynamical range available to 
harmonic coordinates by adding a damping term to the 
equation: 



(A2) 



where is the future directed unit normal to the 
constant-t hypersurfaces. Adding such a damping term 
to the equations for the spatial coordinates tends to 
remove extraneous gauge dynamics and drives the coor- 
dinates toward solutions of the co- variant spatial Laplace 
equation on the timescale Choosing to be com- 
parable to (or smaller than) the characteristic timescale 
of a particular problem should remove any extraneous co- 
ordinate dynamics on timescalcs shorter than the phys- 
ical timescale. The addition of such a damping term in 
the time-coordinate equation is not appropriate however. 



Such a dampcd-wave time coordinate is driven toward a 
constant value, and therefore toward a state in which it 
fails to be a useful time coordinate at all. It makes sense 
then to use the damped-wave gauge condition only for 
the spatial coordinates: 



(A3) 



where N"^ is the shift, and N is the lapse. The appro- 
priate contra-variant version of this damped-wave gauge 
condition is therefore 



Ha = -flsga^N^N, 



(A4) 



where gab = ipab + tati, is the spatial metric.^ 

We point out that the damped-wave gauge condition, 
Eq. (jA4p . is very similar to one version of the F-driver 
shift condition adopted recently by several groups using 
moving puncture evolution methods It is straight- 

forward to express the co- variant wave operator in terms 
of the 3-1-1 decomposition of the metric: 



+g'''dklogN, 



(A5) 



where (^^F* is the trace of the Christoffel connection com- 
puted from gij. It follows that the damped-wave shift 
condition, Eq. ()A3p . is equivalent to the following condi- 
tion on the shift: 



2 (3)-pi 



(A6) 



In comparison a version of the F-driver shift condition, 
cf. Eq. (26) of Ref. that is currently being used by 
a number of numerical relativity groups is a very similar 
condition: 



dtN' - N^dkN' + rjN' = 0.75 ^^^f' 



(A7) 



where (^^F* is the trace of the Christoffel connection com- 
puted from the conformal metric gij = g~^^^gij- This 
version of the F-driver shift condition is therefore a cer- 
tain conformal damped-wave equation for the spatial co- 
ordinates. 

While the damped-wave gauge is a poor choice for the 
time coordinate, the idea of imposing a gauge that uses 
the dissipative properties of the damped-wave equation 
to suppress extraneous gauge dynamics is attractive. The 
lapse is the rate of change of proper time with respect to 
the time coordinate (as measured by an observer moving 
along t"-), so choosing a gauge in which the lapse satis- 
fies a damped-wave equation seems like the appropriate 



^ Frans Pretorius and Matthew Choptuik have recently, indepen- 
dently, proposed adding similar damping terms to the harmonic 
gauge condition Il4ll . 
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time-domain analog of the damped-wave spatial gauge 
condition. To find the appropriate expression for f^Ha 
that leads to such an equation, we note that the gauge 
constraint Ha + Fa = implies that f^Ha is given by 



t'^Ha = -K -f'dalogN, 



(A8) 



where K = g'-' Kij is the trace of the extrinsic curvature 
of the constant-t hypcrsurfaces. Using the definition of 
K , this condition can be also be written in the form, 



= t-da log ( ^ ) - N-^dkN\ 



(A9) 



where g = det gij is the spatial volume element. One fre- 
quent symptom of the failure of simpler gauge conditions 
in binary black-hole simulations is an explosive growth 
in g in the spacctime region near the black-hole horizons. 
This suggests choosing the gauge condition. 



t^Ha = -fiLl0g[^ 



(AlO) 



for /iL > 0, which tends to suppress any growth in ^/g/N 
as a consequence of the constraint, Eq. (jA9p . 

To determine how this gauge condition, Eq. (jAlOp . ef- 
fects the evolution of the lapse, we note that the time 
derivative of K is determined by the Einstein evolution 
equations: 



f'daK 



N-^D'D,N, 



(All) 



where Di is the spatial co-variant derivative compatible 
with gij . Combining this expression with Eq. jASj gives 
an equation for the time derivative of f^Ha, 



Nt'dtit-Ha) 



-t^di,{t''daN) +D'D,N 
+N-\t''daNf - NK,jK'\Al2) 



which is a wave operator acting on the lapse. When the 
gauge condition in Eq. (jAlOp is enforced, it equates this 
wave operator to the following expression, 

Nt''db{eHa) = fiLf'daN - ifiLNt^dalogg. (A13) 



The first term on the right side of Eq. (|A13p is a stan- 
dard damping term for the lapse wave equation, while 
the second term plays the role of an additional "source." 
The motivation for including the particular dependence 
on g in Eq. (|A13|) is provided by the argument leading to 
Eq. (jAlOp . however, this dependence can easily be gen- 
eralized without changing the term's fundamental lapse- 
damping property by setting 



(A14) 



t-Ha = -fiL\0g[2- 



where p is a constant. The case p = 0.5 corresponds to 
Eq. (|A10|) . while p = is a pure lapse-damping gauge 
without the extra source term. 



Combining this new lapse condition, Eq. (|A14p . with 
the damped-wave spatial coordinate condition, Eq. (jA4|l , 
gives the target gauge-source function for our full 
damped-wave gauge condition: 



Fa = ML log ( 2- ) - ^^sN-^ga^N\ (Al5) 



This gauge condition depends only on the spacetime met- 
ric tpabi so it could be implemented directly in the GH 
Einstein system by setting Ha = Fa . However it can also 
be implemented with the new GH Einstein gauge-driver 
system introduced in Sec. |TI1 where it can be used as a 
non-trivial test of the new gauge-driver. Numerical evo- 
lutions of strongly perturbed single black-hole spacetimes 
using the p = 0.5 version of this gauge and the new GH 
Einstein gauge-driver system are described in Sec. IHII 



APPENDIX B: HYPERBOLICITY 

The hyperbolicity of an evolution system consisting 
of some first-order equations, like our new gauge driver 
Eqs. ^ and (HI]), and some second-order equations, like 
the GH Einstein Eq. ([2]), is most easily analyzed by con- 
verting all the equations to first-order form. The spec- 
tral evolution code that we use to perform our numeri- 
cal simulations is rather sensitive to ill-posed evolution 
problems. So we generally perform our numerical simula- 
tions by evolving first-order systems of equations where 
hyperbolicity is easier to analyze and where boundary 
conditions are easier to construct. Mixed systems like 
the combined Einstein and gauge-driver equations can be 
converted to first-order form by introducing additional 
dynamical fields for the first derivatives of those fields 
satisfying second-order equations. Convenient choices of 
the needed additional fields for the GH Einstein system 
are Uab = -fdcipab and = ditpab- The evolution 
equations for these fields, {ipab,^ab,^iab}, then become 
a first-order representation of the GH Einstein system: 

dtyJab - (1 + li)N''dkyJab = -NHab - llN'-P^ab, (Bl) 

dtliab - N'^dkUab + Ng'''^k'^>^ab ' 7l72A^''9fc V-ah 

+2Nd^aHi,) = -^NtH''U,dTlab ~ Nt'n,,g'^^,ab 
+Njo [2S\atb} - ^abt'] [He + r^T,,f) 



(B2) 



^t^^ab - N'^dk^^ab + N d^Xlab - N-f2d^->ljab 

= kNtH''^,,Mab + Ng^h'^^.j.^kab - Af72*«fc, (B3) 

cf. Eqs. (35)-(37) of Rcf. 0. In these equations TV, N\ 
and gij, are the standard 3-1-1 representation ofipab given 
in Eq. ([7]); t° is the future directed timelike unit normal; 
r^j, is the Christoffel connection associated with Tpab] and 
7o, 7i, and 72 are parameters multiplying constraints, in- 
troduced because they help damp away small constraint 
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violations. This representation of the GH Einstein equa- 
tions together with the gauge driver introduced above, 
Eqs. (O and pTjl . is a first-order evolution system which 
can be represented abstractly as, 



The characteristic speeds. 



•J{a): 



associated with the 



(B4) 



For the combined GH Einstein gauge-driver sys- 
tem, the collection of dynamical fields is = 
{ipab, Hab, ^iab. Ha, 0a}, where Greek letters are used for 
indices that enumerate the dynamical fields. 

The hyperbolicity of a first-order evolution system, 
such as Eq. (|B4p . is determined by the properties of the 
characteristic matrix A^^-p. We define the left eigen- 
vectors e"^ and their associated eigenvalues V(^a) of the 
characteristic matrix in the following way, 



(B5) 



where rife denotes a spacelike unit vector; accented Greek 
letters, d, are used to enumerate distinct linearly in- 
dependent eigenvectors. The eigenvalues, W(q,), are often 
referred to as the characteristic speeds of the system. A 
first-order evolution system is strongly hyperbolic at a 
point in spacetime if there exists a complete set of eigen- 
vectors for each Uk at that point. In this case the matrix 
of eigenvector components e"a is non-degenerate, i.e., 
det e^Q ^ 0. The projections of the dynamical fields onto 
the eigenvectors, = c^qm", provide an alternate com- 
plete set of dynamical fields, which play an important 
role in strongly hyperbolic systems. For example, the 
characteristic fields, u", are those on which appropriate 
boundary conditions must be placed for these systems. 

It is fairly straightforward to work out the characteris- 
tic eigenvalues and eigenvectors, and the associated char- 
acteristic fields, for the combined GH Einstein gauge- 
driver system: 



where P, 



^ab, (B6) 
Uab ± n^^iab - J2lpab ± naHb ± UbHa, (B7) 
P^'^^kab, (B8) 
Ha, (B9) 
ea+VHa, (BIO) 

Si^—riiv}' . We see that the coupling between 
the GH Einstein and gauge-driver systems increases the 
number of characteristic fields, and also transforms the 
characteristic fields of the pure GH Einstein system. This 
means that the theory of the boundary conditions for 
the GH Einstein system will have to be completely re- 
examined. We also note that the co-vector Ua is a spa- 
tial unit normal, which is orthogonal to the timelike unit 
normal ta- This implies that the spatial components of 
Ua are the usual components of the spatial normal co- 
vector Ui while the time component rit must be given by: 
rit = fikN'^. These conditions ensure that t°n„ =; and 



combined GH Einstein gauge-driver system are as fol- 
lows: the fields have coordinate characteristic speed 
-(l+7i)7ifeiV'=, the fields have speed -nkN^±N, the 
fields M^^j and uj^ have speed —ntN^, and the fields u\ 
have speed zero. On boundary points each characteristic 
field (computed with the outward directed unit normal 
to the boundary rik) must be supplied with a boundary 
condition if and only if its associated characteristic speed 
is negative. The appropriate boundary conditions for the 
combined GH Einstein gauge-driver system are discussed 
in some detail in Appendix [El 

The inverse transformation between dynamical and 
characteristic fields for the combined GH Einstein gauge- 
driver system is 



V'oh 
Hah 



-ni{naul + Ubul), 



Ha = K 



(Bll) 
(B12) 

(B13) 
(B14) 
(B15) 



Since this transformation is invertible, the combined 
first-order GH Einstein gauge-driver evolution system is 
strongly hyperbolic. 



A first-order evolution system, Eq. (|B4[) , is called sym- 
metric hyperbolic, if there exists a symmetric positive 
definite matrix on the space of dynamical fields, Sais, that 
symmetrizes the characteristic matrices: Sa-yA^'^ p = 
A'^P = Ap^. Symmetric hyperbolic systems provide 
a natural "energy," E = J Sapu^u^d^x, and arc bet- 
ter behaved than strongly hyperbolic systems for initial- 
boundary value problems. Symmetric hyperbolicity is 
therefore a desirable property for gauge-driver systems 
to have. It is fairly straightforward to show that the 
combined GH Einstein gauge-driver system of Eqs. ([2]), 
^ and pT|) has a symmetrizer given by: 



dS'^ = S^pdu^duf^, 



A^dHadHb 



I ^o,b^cd r ij 

+m m [g 



AlidOa + vdHa){d0b + vdHb)] 
{d<^^ac + "igiadHc) {d^jbd + 2gjbdHd) 

{dliac - ^2dlpac) (dUbd - Jldlpbd)] .(B16) 



1. 



This symmetrizer is positive definite as long as to"'' is a 
positive definite symmetric tensor, and the (real) scalars 
A^, Ah, and Ag are non- vanishing. Therefore the gauge- 
driver system of Eqs. © and ITTI) is symmetric hy- 
perbolic. 
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APPENDIX C: DUAL COORDINATE FRAMES 

We have found that using two different coordinate sys- 
tems simultaneously is a very useful numerical technique, 
when performin g nu merical evolutions of binary black- 
hole spacctimes [IJ]. This method allows us to choose 
one set of coordinates, thought of as "co-moving," to 
track (approximately) the motion of the black holes, and 
a second set, a;" thought of as "inertial," fixed (approx- 
imately) to a non-rotating frame at infinity. We evalu- 
ate the components of the various dynamical fields using 
tensor bases defined by the inertial x"" coordinates, while 
the evolution equations arc solved for those inertial-frame 
field components as functions of the moving x'" coordi- 
nates. This use of dual coordinate frames minimizes the 
size of the various field components and their time deriva- 
tives better than any single-frame coordinate choice. 

The single- frame GH Einstein gauge-driver equations, 
introduced in Sec. HIl written in terms of inertial-frame 
quantities are given by 



-HiH-a - F-a) + e-a, (CI) 
-VO-a- (C2) 



These equations, together with the inertial-frame repre- 
sentations of the Einstein system, can be converted to 
dual-frame form in a straightforward way using the pre- 
scription developed in Ref. [lj|. Under this recipe, a first- 
order evolution system for inertial frame components, m". 



is converted into the dual-frame system 



5m" 



k a _ 



(C3) 



(C4) 



simply by changing independent variables: dt = dt + 
djx^di and dj, = dj.x'^di. The quantities dix^ = dx^ jdt 
and d^x"^ = dx^ jdx^ are the non-trivial parts of the Ja- 
cobian of the transformation relating the two coordinate 
frames. These coordinate transformations are assumed 
to be given a 'priori. 

The straightforward conversion of the GH Einstein 
gauge-driver system from its inertial single-frame form, 
(jCip . and (|C2[) . to dual- frame form may not always be 
the most effective choice however. The single-frame evo- 
lution equation for i/g, Eq. (|C1|) . is designed to drive 
Ra Fa at fixed values of the inertial coordinates. 
A binary black-hole spacetime, however, can have rapid 
time variations in the field components when evaluated 
at fixed inertial coordinates, e.g., at points lying near 
the black-hole trajectories. The gauge-driver system will 
not be very efficient in accurately enforcing the desired 
gauge under these very dynamical conditions. In contrast 
the moving coordinates, x"", are chosen to track (approx- 
imately) the motion of the holes, so the fields expressed 
as functions of the moving coordinates are far less time 
dependent. A moving- frame version of the gauge-driver 



would therefore be more effective enforcing the desired 
gauge. Ha = Fa, in many situations. In this case it 
makes sense to modify the evolution equation for Ha in 
a way that ensures the moving-frame components of Ha 
are driven to the intended targets: Ha — > Fa- The ap- 
propriate moving-frame gauge driver equations are sim- 
ply Eqs. (jlip and ^ interpreted now as moving-frame 
equations: 



dtHa - N'dkHa = -fi{Ha - Fa) 
dtOa+vN^duHa^-llOa. 



(C5) 

(C6) 



It is straightforward to re-express these equations in 
terms of inertial frame quantities: 



dtH-a - N'^d-^H-a = -^liH-a - Fa) + 9a 
+ (d-tdaX'' ~ N^'d-^d-aX'^) dax'Hj,, 



(C7) 



diO-a + dtX^d^e-a + 77 (iV^ + dtx'') d-kH-a = -TJ Oa 
+ (af5ax" + dtx''dj,daX'')dax''0-^ 

+Tj{N^ + dtx'^){dj,daX'')dax''H-„ (C8) 

where da = daX°'da transforms the derivatives. Ha = 
dax"'Ha and 9a = daX°'9a transform the field components, 
and the inertial-frame shift N'^ is related to the moving- 
frame shift N'' by 



(C9) 



In some circumstances it may be advantageous to 
apply the inertial-frame version of the gauge driver, 
Eqs. (jCip and (jC2p . in one region of spacetime, while 
applying the moving-frame version, Eq. (|C7|) and (|C8|) . 
in another. For example, in a binary black-hole simula- 
tion it might be appropriate to impose the moving-frame 
version of the gauge driver in the very dynamical region 
of spacetime near the black holes, while imposing the 
inertial-frame version in the more quiescent asymptotic 
region far from the holes. Therefore, to accommodate 
this possibility we introduce the following hybrid gauge- 
driver system that simply interpolates between the two: 



diH-a - N^d-kH-a = -f4Ha - Fa) + Oa 

+W (d-td-aX"- - N'^d-^d-aX'') dax'H-t,, 



di9a + wdtx''d-k9a + r/ [N^ + wdtx^j d-kH-a 

^widid-ax" ^ dtx^d-kd-ax'')dax^9i 
+wr,{N'' + dtx'')idkdax'')dax''H-^. 



(ClO) 

-VO-a 

(Cll) 



In these equations the smooth weight function w is spec- 
ified a priori, with w = in the spacetime region where 
an inertial-frame gauge driver is needed, and w = 1 in the 
regions where a moving-frame gauge driver is required. 
The dual-frame version of this hybrid gauge-driver sys- 
tem is obtained by combining these equations with the 
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inertial- frame Einstein system, Eqs. (|Bip - (|B3p . and us- 
ing the dual-frame conversion technique summarized in 
Eqs. jCH and (IC4ll . 

Since the hybrid gauge driver Eqs. (jClOP and (jClip 
do not have the same principal parts as their single- 
frame counterparts, we must consider again the hyper- 
bolicity of the combined GH Einstein plus hybrid gauge- 
driver system. Fortunately, wc find that this system is 
still strongly hyperbolic, and the characteristic fields are 
just Eqs. (|B6|1 - (|B10|1 expressed in terms of inertial-framc 
field components. The characteristic speeds associated 
with these fields are modified somewhat however: The 
fields u^j^ have inertial-coordinate characteristic speed 

-(1 + 7i)rijiV^, the fields u^f have speeds -n^N^ ± N, 

the fields and have speed —nj.N'', and has 

the speed wrif.dtx''. In these expressions N and N'^ refer 
to the inertial frame lapse and shift respectively. The 
co-moving-frame characteristic speeds are obtained from 
the inertial-frame speeds by adding ~n^.dtx^ . This hy- 
brid gauge driver system is also symmetric hyperbolic 
with the same symmetrizer, Eq. (|B16p . interpreted as an 
expression in terms of inertial-frame field components. 



APPENDIX D: CONSTRAINTS 

This appendix investigates the constraints of the new 
GH Einstein gauge-driver system. These constraints and 
(somewhat surprisingly) their evolution equations turn 
out to be identical to those of the pure GH Einstein sys- 
tem. This means that the constraint-preserving bound- 
ary conditions derived for the pure GH Einstein system 
are also appropriate for the combined GH Einstein gauge- 
driver system, although care must be taken to enforce 
them on the correct characteristic fields of the combined 
system. This section presents the groundwork for the de- 
tailed discussion of boundary conditions in Appendix lEl 

The primary constraint of the GH Einstein system is 
the gauge constraint, Ca, which can be written in terms 
of the first-order dynamical fields: 

(Dl) 

There are no extra constraints from the addition of the 
first-order gauge-driver fields Ha and 9a- In the pure GH 
Einstein system the gauge-source function Ha is assumed 
to be a prescribed function of the spacetime coordinates 
x°' and the 4-metric ipab'- Ha = Ha{x,ip). In contrast 
Ha is elevated to the status of an independent dynamical 
field that is evolved according to Eq. (fTTjl in the combined 
GH Einstein gauge-driver system. We need to determine, 
whether the evolution of the GH constraint fields is af- 
fected by the introduction of this gauge-driver equation. 
In addition we need to find the characteristic constraint 
fields to determine what constraint preserving boundary 
conditions are needed for the new combined system. 



The basic GH Einstein system, Eqs. (jBip (jB3p is, as 
before, just a representation of the 4-dimcnsional co- 
variant Einstein equation: 

Rab = ^iaCb) - 70 [haCb) ^ i^Jabt'Cc] , (D2) 

where Rab is the Ricci curvature, and Va is the co- variant 
derivative associated with il'ab- Consequently the evolu- 
tion equation for Ca is determined by the Bianchi iden- 
tities for the 4-dimensional Ricci tensor, which can be 
written as the second-order wave equation: 

= WcCa~2-/oy%bCa)]+C''y^aCb) - ilotaC'Cb. 

(D3) 

This equation is identical to that obtained for the pure 
GH Einstein system 0], because its derivation does not 
depend on how the Ha field is evolved. 

Constraint preserving boundary conditions are de- 
signed to prohibit the influx of constraint violations 
through the boundaries of the computational domain. In 
order to fix the incoming constraint fields, the charac- 
teristic fields of the constraint evolution system must be 
identified. This is done most easily by transforming the 
second-order constraint evolution Eq. (jPSp to first-order 
form. To do this we introduce new constraint fields rep- 
resenting the first derivatives of Ca- Thus we define new 
constraint fields J-'a and Cia that satisfy 

Ta « t'dcCa = N-\dtCa - N'O^Ca), (D4) 
C^a ~ d,Ca, (D5) 

where « indicates equality up to terms proportional to 
the gauge constraint Ca and the first-order GH Einstein 
constraint Ciab = SiV'ah ~ ^iab- The following expres- 
sions for J^a and Cia accomplish this in a way that keeps 
the form of the constraint evolution system as simple as 
possible: 

Ta = i5^^^^a,H,,-5^^a,H,,-5*^t^a,;$,fc, 

+gi^ijhg'''^kcd^'"'t' - igl^^jbd'^^kcdr't" 

+ ^ta^'''9''d,^jbc - \tag''<^^cd'i>jbe^'V 
-ita.g''^5""$,™c$„,dV'''' + g'''^>^cd<^Jbai^''t^ 
+ itaHednf,eV''V''' " 9^' H.Ilja ~ t^'g^mullja 
-ig^^.dt^t'^HbeV''' + itaU^dTlteij^HH' 
+g'a^^cdllbetH''^''' - ^g'^ <P,,dtH''U,a 

-g'^^^bat'Tl^et" + gl<S>^cdHb^"=t'' + tag'^d.H, 
+i2{g'%da - igli^^'^C.cd) + ktag^'H^^.^dr" 

+ ^,taIlcdi'"'Hbt'' - g'^H.^jbat' - glt'd.Hb 

-tag'^'^,j,Hdr\ (D6) 

C^a = g'''dj^^ka - kgi^^'^dj'^^cd - kta^P'''^^n,d 
+t''dinba + d,Ha + hgi'^jcd'^^ef^'^V^ 
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+ hl2 {tar" - '^Slt") C,cd. (D7) 

We note that while Ta is defined as the time deriva- 
tive of Cq, the expression in Eq. (|D6p contains no time 
derivatives. The constraint fields are functions of the fun- 
damental dynamical fields of the system m". Any time 
derivatives of the constraint fields are determined by the 
time derivatives of these fundamental fields through the 
evolution equations of the system. When the time deriva- 
tives of the expression for Ca in Eq. (|Dip are evaluated, 
and the time derivatives of {ipab^^ab, '^iab\ are replaced 
with the expressions from the basic GH Einstein system, 
Eqs. (jBip - (|B3p . we find that the occurrences of dtHa 
cancel one another. Thus the expression for Ta does not 
depend on how Ha is evolved, and it is valid for both the 
pure GH Einstein system and the new first-order gauge 
driver system. To complete the GH constraint evolution 
system we need to add the GH Einstein constraint Ciab^ 

C^ab ^ d,l}jab - ^tab, (D8) 

and the closely related Cijab, defined by 

Cijab = 29[j$j](j(, = 2d[jCi]ab- (D9) 

The complete collection of constraints for the GH 
Einstein gauge-driver evolution system is the set = 

{Ca,J'a,C^a,C^ab,C^Jab} defined in EqS. (|D7|, 

(|D8p . and (|D9|) . (We use upper case Latin indices to enu- 
merate the constraint fields.) The constraints depend 
on the dynamical fields = {ipab, ^ab, ^iab, Ha, 9a} and 
their spatial derivatives dku". We have evaluated these 
constraint evolution equations using the new GH Einstein 
gauge-driver system and have verified that they can be 
written in the abstract form 

dtc' + A'' ' j{u)dkc'' = j{u,du) c\ (DlO) 

where A''^ j and j may depend on the dynamical fields 
and their spatial derivatives dku°'- The evolution of 
the constraint fields turns out to be completely deter- 
mined by the GH Einstein Eqs. (jBip (|B3p alone without 
any use of the gauge-driver Eqs. ([9]) and (|lip . While 
the constraint fields Ca, Ta and Cia depend on Ha and 
dkHa, the time derivatives of these constraints are de- 
termined without using the evolution equation for Ha, 
Eq. pT|) . There is a remarkable cancellation between the 
explicit time derivatives of Ha appearing in dtTa and 
dtCia, and the time derivatives of Ha introduced when 
the dfllta terms arc replaced in these expressions using 
the GH evolution Eq. (|B2p . Thus the constraint evolu- 
tion system for the first-order gauge-driver system does 
not depend at all on the gauge driver Eqs. ([9]) and (fTTj) . 
This constraint evolution system is identical to the pure 
GH Einstein constraint evolution system given in Ref . , 
and is both strongly and symmetric hyperbolic. 

Since the constraint evolution equations for the GH 
Einstein gauge-driver system are identical to those of the 



pure GH Einstein system, the characteristic constraint 
fields are also identical. The boundary conditions 
needed to ensure no influx of constraint violations will 
also be the same therefore. As we have seen in Ap- 
pendix [B] however, the characteristic dynamical fields u" 
of the two systems are not the same, so the detailed 
expressions for the needed boundary conditions in the 
two systems will be different. So we recall here the ex- 
pressions for the characteristic constraint fields from 
Ref. 0: 

4* - :FaTn''Cka, (Dll) 
cl ^ Ca, (D12) 
4 = P^Cka, (D13) 

cfab = C,ab, (D14) 
(^ijab ^ Cijab- (D15) 

The characteristic constraint fields have coordinate 
characteristic speeds —niN''ztN, the fields have speed 
0, the fields and cfjab have speed —niN'-, and the fields 

'^iab have speed ~{l+ji)niN'' . Boundary conditions must 
be placed on the incoming characteristic dynamical fields 
that (among other things) fix the incoming charac- 
teristic constraint fields to zero. These (and other) 
needed boundary conditions are discussed next in Ap- 
pendix. [e1 

APPENDIX E: BOUNDARY CONDITIONS 

A boundary condition is required for each characteris- 
tic field of the GH Einstein gauge-driver system, at 
each boundary point where the characteristic speed W(q,) 
associated with that field is negative. 

The characteristic fields Eq. (jB6|) . have speed 

— (1 -|- "fi)nkN^ and may require boundary conditions 
at some boundary points. Since the constraints and 
the constraint evolution equations of the GH Einstein 
gauge-driver system are identical to those of the pure GH 
Einstein system, we can employ the same approach to 
constructing constraint preserving boundary conditions. 
The constraint characteristic field c^^^, Eq. (jD14p . is re- 
lated to the characteristic field vPab hy the expression, 

n'cl, « d^ul,, (El) 

where dj^u" denotes the characteristic projection of the 
normal derivatives of u", i.e., d±u°' = pn'' dkU^ , with 
e"/3 defined in Eq. (jBSp . Here (and throughout this ap- 
pendix) w implies that algebraic terms and terms involv- 
ing tangential derivatives of the fields (e.g. Pi^dkU°') 
have not been displayed. We note that the constraint 
field cfab has the same characteristic speed as Hence 
a constraint preserving boundary condition for c^^^ is 
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needed whenever needs a boundary condition. The 
identity relating to cf^^, Eq. (|Eip . provides the way 
to formulate this boundary condition by prescribing the 
value of d_LW°j,. 

A convenient way has been found Q to impose con- 
straint preserving boundary conditions for fields like uj^^ 
that are related to an incoming constraint field through 
an expression like Eq. (jEip . The characteristic projection 
of the time derivatives of these fields u", dtu" = e"pdtu'^, 
are set in the following way at the boundary, 

dtu" = Au" + U(a)(d_LM"-rf±u"|BC')- 

In this expression the terms Dfu" represent the projec- 
tions of the right sides of the evolution system, Eqs. (|B1[) - 
(|B3p . (|ni) and pi|l . so the characteristic projections of the 
evolution equations at non-boundary points would sim- 
ply be dtw" — Dfu". The term d±u"\g^ is the value to 
which d±u°' is to be fixed on the boundary. This form 
of the boundary condition replaces all of the d^u" that 
appears in Dfu" with d±u"\^^. Applying this method 

to the u^fj field, wc arrive at the desired constraint pre- 
serving boundary condition for this field, 

dtut, = Dtui, - (1 + 7i)njiV^n^cL. (E3) 

This boundary condition is the same in the new GH Ein- 
stein gauge-driver system as in the pure GH Einstein 
system 0. 

The characteristic field uf^f,, Eq. (jBSp . has speed 
— nfeTV'^, and so this field may require a boundary con- 
dition on some boundary points. The constraint charac- 
teristic field c^jjjft, Eq. (jPlsp . has the same characteristic 
speed, and hence it is natural to use the boundary con- 
dition on uf^fj to prevent the infiux of this constraint. 
Conveniently, there is an identity relating u'^^^^ and cfj^f^: 

n'cjkab ~ d^ul^^. (E4) 

This identity is identical in the GH Einstein gauge-driver 
and the pure GH Einstein systems [?[■ So we follow the 
strategy of Ec^. (jE2p . and use the following constraint 
preserving boundary condition for wf^f,: 

dtut, = Dt4ab-niN'n'Ph4ab- (E5) 

The characteristic field u^, Eq. (|B9p . also has speed 
—UkN'^, and so it may require a boundary condition on 
some boundary points. We have identified two possibil- 
ities for this boundary condition. First, the Ha field is 
part of the basic gauge constraint Eq. (jPip . So one pos- 
sible boundary condition for is simply to enforce this 
constraint on the boundary: 

"a = -.9''*ya - t'Hba + k_ gl^''' <i> ^bc + ^iaV^'^Hfee- 

(E6) 



Another possibility is to use a boundary condition on 
that enforces the desired gauge condition Ha = Fa on 
the boundary: 

4 = Fa- (E7) 

These boundary conditions could be imposed as Dirichlet 
conditions in the forms given above using penalty meth- 
ods. Alternatively, we could impose these conditions us- 
ing Bjorhus methods as a driver condition on the bound- 
ary value of the time derivative of the characteristic field, 

dtu" = -Mu" -u"\bc)- (E8) 

The parameter fj,B sets the timescale on which the bound- 
ary value of u" is driven to its target value. The Bjorhus 
version of the boundary condition in Eq. (|E6P is there- 
fore, 

dtul = ~^iBCa, (E9) 

while the Bjorhus version of Eq. (|E7p is 

dtul = -tiB{Ha-Fa). (ElO) 

In most of our numerical tests, we find that the Eq. (jElOp 
version of this boundary condition is more effective. 

The characteristic field m^, Eq. (jBlOp . has character- 
istic speed in the single-frame evolution system, and 
hence does not need a boundary condition in that case. 
In the dual-frame system the characteristic speed changes 
to wnj^dtx'^, so this field might need a boundary condi- 
tion under some conditions. We have generally chosen 
weight functions w and dual-frame maps dtx'' that avoid 
the need for a boundary condition on this field. But if 
that can not be done, it is probably best to choose the 
boundary value of 6a so that dfOa = on the bound- 
ary. This condition leads to the following Dirichlet type 
boundary condition for u^: 

ul^f^Ha-N''dkHa. (Ell) 

Boundary conditions are rarely needed for the 
fields, i.e., only when the boundary of the computational 
domain moves outward at supcrluminal speeds. In con- 
trast boundary conditions arc almost always needed for 
the it^j^ fields. These boundary conditions split naturally 
into three types that have been called gauge boundary 
conditions, constraint-preserving boundary conditions, 
and physical boundary conditions 0, These three 

different types of boundary conditions are imposed on 

the parts of selected by the three mutually orthogo- 
nal projection tensors: 

Pab^"' = -[kahl'^' + kaSb^' + hSj']!''^ (E12) 
PaP""" = iPa6P'=''-2Z(„Pfc)(^fc'^)+U6A:^fc^(E13) 
^ir^^ = Pj'Pb''^ - iPabP^''. (E14) 
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In these expressions fc° and Z° represent the ingoing and 
outgoing null vectors, respectively, that are related to the 
timelikc and outgoing spacelike unit normal vectors, 
and by: 



^/2 V 



(E15) 
(E16) 



Similarly gat represents the spatial 3-metric and Pab the 
projection onto the 2-dimensional spatial boundary sur- 
face: 



gab = 'ikab + tah, 

Pab = i^ab + tah- riaUb. 



(E17) 
(E18) 

Finally, we note that the projection tensors Pj^^'''^'^, 
P^b^'^'^^ and P^b^'^'^ are complete in the sense that: 

Sa^'Sb'''> = P^^f^'"" + P^^^'^ + Pi^""- (E19) 

We now discuss the boundary conditions appropriate for 
the three independent projections of the u]^ fields. 



So using Eq. (jE2ip we set: 



\BC 



= {Sa' - tat') {-Hb - g''<^,,b + igbH'^'^.cd 



+ it6.g'-'n.y) 
Uabt' - {Sa" - tat')Cb. 



(E22) 



Using this expression in the equation for u^^^ l'' 
Eq. (|E20p , we find the expression for the target boundary 
value of u^^l'' to be: 

^bl'lsc = <^l'-^{S,'-tat')Cb. (E23) 

This boundary condition can either be imposed as 
a Dirichlet condition by penalty methods, or as a 
boundary-driver condition by Bjorhus methods using 
Eq. (jESp . The Bjorhus version of this new gauge bound- 
ary condition is: 



p(G)cd 7 ,,i 
^ab 



t^d = -^'BPab 
MB 



(G)crf/ i- 



V2 



kn.kb^ 



kaSb 



'■cd \bCI 
+ kb5n 



X {5/ - tct'')Cd. 



(E24) 



2. Constraint Preserving Boundary Conditions 



1. Gauge Boundary Conditions 



The term gauge boundary conditions is used to de- 
scribe the boundary conditions on the P^^ projection 



of ^ 



From the structure of the P, 



iG}cd 



projection 



tensor, we see that these are in effect boundary condi- 
tions on the u^^^l'' fields. Writing out the definition of 



^ab ' that 



ab 



Ilabr - n'^^abl' - ^2la ' TlaHbr 



(E20) 



The u]^ characteristic fields determine the time and the 
spatial derivatives of ■0a6 normal to the boundary. So 
these gauge boundary conditions on u]~^l^ can be thought 
of as fixing the liab^' components of Ilafc. Previously 
the gauge boundary condition on u]~^l'' has been set by 
freezing the value of this projection of the characteristic 
field, Pah^'"^dtu\'^ = 0, or by imposing a Sommerfeld- 

like boundary condition on this projection of u]~^ WA - 

Here we present a new gauge boundary condition for 
u\'^l' obtained by setting the target boundary value of 
Ti-abt' to the value it would have if the gauge constraint 
were satisfied exactly. The components of Ilabt' enter the 
gauge constraint, Ca, through the identity: 

nafct' = {5a'' - tat') {Cb ^ Hb - .g'^ ^yb + i^&W^^'^.cd 

+itbg'm,j). (E21) 



The term constraint preserving boundary conditions is 
used to describe the boundary conditions on the P^'j^^'^'^ 

projection of u^b ■ These boundary conditions have 
been constructed to enforce the incoming components 
of the constraint characteristic fields c° = 0, defined in 
Eq. (jDll|) . at the boundary. For the pure GH Einstein 
system, it was shown that 



cdl 



xn'^dkiUcd 



n'^'icd -7oV'cd).(E25) 



For the case of the pure GH Einstein system, this gives 
an expression for c^^ in terms of the normal derivative 
and so can be used to construct a boundary con- 



of 1 

dition using Eq. (|E2|) . For the new GH Einstein gauge 

driver considered here, the u^^ characteristic fields in- 
clude the additional terms —ii(^aHb). In the derivation 
of Eq. (|E25)) from Eqs. ((D6)) and (|D7)) . the terms in- 
volving spatial derivatives of Ha were treated as being 
prescribed, and so were counted as some of the (many) 
algebraic terms not displayed. Since Ha has been ele- 
vated to the status of a dynamical field in the new first- 
order gauge-driver system, however, these terms can no 
longer be ignored. It turns out that the dkHa terms in 
Eqs. (jD6p and (|D7p give the following extra contributions 
to Eq. (|E25l) : 



V2[k^'i/\ - ^^ka^' 



UrHd - TldHc) 

(E26) 
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Using this expression in Eq. (|E2p . we then arrive at the 
needed boundary condition for the constraint preserving 
components of : 

X [haPb)' - ^Pabl' - i^a/^,fc"]c^. (E27) 

These boundary conditions have the same form as those 
derived previously for the pure GH Einstein system Q. 

Here however, the characteristic field u]^'^ has a different 
meaning, since it depends explicitly on the Ha field in 
the GH Einstein gauge-driver case. 

3. Physical Boundary Conditions 

The term physical boundary condition is used to de- 

(p\ , 

scribe the boundary condition on the P^j, projection 

of 01 ■ This projection corresponds to the transverse 
traceless components of the metric field, and so describes 
the physical gravitational wave degrees of freedom of the 
system. In the vacuum region far away from compact 
sources, the gravitational- wave degrees of freedom are de- 
scribed by the propagating components of the Weyl cur- 
vature tensor. The characteristic fields, w^^^, representing 
these incoming and outgoing wave degrees-of-freedom re- 
spectively of the Weyl tensor, arc given by 

<b = P!,P'''{t'Tn^WTnf)C,edf. (E28) 

It is straightforward to show that the incoming 
gravitational-wave characteristic field w"^ depends on the 
normal derivatives of the dynamical fields at the bound- 
ary by the expression: 

Kb « PiP'^n'^d,{n,d-n^<^,,,i), (E29) 

where w denotes that algebraic terms and terms depend- 
ing on tangential derivatives of the dynamical fields are 
not shown. The derivation of this expression depends 
on the fact that the physical projection P^^^ ' annihi- 
lates terms like PiP'^^i^cd = and P^P'^'^n^cHd) = 0. 



Therefore the principal part of w^i, depends on the nor- 
mal derivative of u^]^: 

Kb « PiP^'d^^d- (E30) 

This is the same expression (up to terms proportional 
to constraints) that is satisfied in the pure GH Einstein 
system case [7| , where the gauge-source functions are pre- 
scribed: Ha = Ha{x, But here the characteristic field 
has a somewhat different meaning. The lowest-order 
physical boundary condition is designed to enforce the 
no-incoming wave condition u)^^ = at the boundary. 
It docs this by using Eq. (jE30p to replace the normal 
derivative of u^j^ which appears in the Einstein evolution 
equation for . This boundary condition is enforced as 

a Bjorhus condition on , 

PtP^'d^ul- = P!,P^'iD,ul- -{N + n,N'')w;,] , (E31) 

which is the same condition used in the pure GH Einstein 
system case 0. Higher-order physical boundary condi- 
tions have also been derived for the pure GH Einstein 
system [l^ , and these could be used, essentially without 
modification, for the GH Einstein gauge-driver system as 
well. 
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